{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Eignenvalues and eignevectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definitions\n",
    "\n",
    "Given a $n\\times n$ matrix $A$, the eigenvalue problem consists in finding a scalar $\\lambda$ and a non-trivial vector $x$ such that \n",
    "\\begin{equation}\n",
    "(E):~~Ax = \\lambda x.\n",
    "\\end{equation}\n",
    "We say that $\\lambda$ is an eigenvalue of $A$ and $x$ is its associated eigenvector. \n",
    "\n",
    "Notice that $x$ is not unique since $\\alpha x$ with $\\alpha\\neq 0$, is still an eigenvector for $\\lambda$. Moreover, if $\\lambda$ s known, the associated eigenvector can be recovered using the so-called *Rayleigh quotion*\n",
    "$$\n",
    "\\frac{\\bar{x}^{T}Ax}{\\Vert x\\Vert^2}\n",
    "$$\n",
    "\n",
    "Studying eigenvalues and eigenvector is of the main interests in many problems in applied mathematics since it allows us to understand the behaviour of the linear transformation $A$. Indeed, equation (E) states that $A$ has a scaling effect on $x$, i.e., flipping, streching or compressing.\n",
    "\n",
    "From now on, assume that $A$ has real entries with real eigenvalues.\n",
    "### Plotting vectors\n",
    "\n",
    "We define a function that takes as inputs a vector $x$ and a matrix $A$ and plots the vectors $x$ and $Ax$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Complete this code\n",
    "def plot_vector(x, A, xlim, ylim):\n",
    "    \"\"\"\n",
    "    function to plot two vectors,\n",
    "    x - the original vector\n",
    "    A - the linear transformation\n",
    "    xlim - the limit for x\n",
    "    ylim - the limit for y\n",
    "    \"\"\"\n",
    "    plt.figure(figsize = (10, 6))\n",
    "    #.....\n",
    "    plt.xlim(xlim)\n",
    "    plt.ylim(ylim)\n",
    "    plt.xlabel(\"X\")\n",
    "    plt.ylabel(\"Y\")\n",
    "    plt.legend()\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your function with\n",
    "\\begin{equation}\n",
    "A = \\begin{pmatrix} 2& 0 \\\\ 0 &1\\end{pmatrix},~x = \\begin{pmatrix} 1\\\\ 1\\end{pmatrix}.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Power method\n",
    "\n",
    "Consider $A\\in\\mathbb{R}^{n\\times n}$ and assume that it has $n$ real eigenvalues $\\lambda_1,\\dotsc,\\lambda_n$ with corresponding eigenvectors $v_1,\\dotsc,v_n$ which are assumed to be linearly independent. Suppose that\n",
    "$$\n",
    "\\vert\\lambda_1\\vert > \\vert\\lambda_2\\vert\\geq \\dotsc\\geq \\vert\\lambda_n\\vert\n",
    "$$\n",
    "and pick some initial vector $x_0$. Since the family $(\\lambda_i)_{i=1}^{n}$ is linearly independent, we can write\n",
    "$$\n",
    "x_0 = \\sum_{i=1}^{n} a_i v_i~\\mbox{with}~a_1\\neq 0.\n",
    "$$\n",
    "Then, multiplying by $A$ in both sides we get \n",
    "$$\n",
    "Ax_0 = \\sum_{i=1}^{n} a_i Av_i = \\sum_{i=1}^{n} a_i \\lambda_i Av_i  = a_1 \\lambda_1 \\underbrace{\\Big(v_1 + \\sum_{i=2}^{n}\\frac{a_i\\lambda_i}{c_1\\lambda_1}v_i\\Big)}_{x_1} = a_1\\lambda_1 x_1\n",
    "$$\n",
    "Again, applying $A$ to $x_1$, we get\n",
    "$$\n",
    "Ax_1 = \\lambda_1v_1 + \\sum_{i=2}^{n}\\frac{a_{i}\\lambda_{i}^{2}}{c_1\\lambda_1}v_i = \\lambda_1\\underbrace{\\Big(v_1+\\sum_{i=2}^{n}\\frac{a_{i}\\lambda_{i}^{2}}{c_1\\lambda_{1}^{2}}v_i\\Big)}_{x_2}\n",
    "$$\n",
    "Continuing this process, we obtain at iteration $k$\n",
    "$$\n",
    "Ax_{k-1} = \\lambda_1\\underbrace{\\Big(v_1+\\sum_{i=2}^{n}\\frac{a_{i}\\lambda_{i}^{k}}{c_1\\lambda_{1}^{k}}v_i\\Big)}_{x_k}\n",
    "$$\n",
    "Since $\\lambda_1$ is the dominant eigenvalue, we have that $\\lambda_i/\\lambda_1 < 1$ for all $i>1$, and consequently $(\\lambda_n/\\lambda_1)^{k}$ goes to zero as $k\\to \\infty$. That is \n",
    "$$\n",
    "Ax_{k-1}\\sim \\lambda_1 v_1.\n",
    "$$\n",
    "Usually, when implementing this method, the fector obtained at each iteration is normalized.\n",
    "\n",
    "> **Exercise**: Implement the power method at test it with\n",
    "$$\n",
    "A = \\begin{pmatrix} 0& 2 \\\\ 2 &3\\end{pmatrix},~\\mbox{and}~A = \\begin{pmatrix} 1& 1 &0 \\\\ 1 &1&1\\\\0&1&1\\end{pmatrix}.\n",
    "$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Insert your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare your result to the function `scipy.linalg.eig`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Exercise**: Assume that $A$ is invertible. Exploiting the implemented power method, compute the dominant eigenvalue of $A^{-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Exercise**: What are the possible techniques to speedup the computations in the previous question $A^{-1}$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:monenv]",
   "language": "python",
   "name": "conda-env-monenv-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}